pkgs <- c("tidyverse", "tidylog", "sp", "raster", "devtools", "RCurl", "sdmTMB", "terra", "ncdf4", "chron")
if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Source code for map plots
source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
# Packages not on CRAN
# devtools::install_github("seananderson/ggsidekick") # not on CRAN
library(ggsidekick)
theme_set(theme_sleek())
# Set path
home <- here::here()Make the prediction grid
Intro
Make an evenly spaced UTM prediction grid with all spatially varying covariates for the diet and the biomass data
Read data and depth-raster
# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") |>
rename(X = x, Y = y)Rows: 9376 Columns: 26
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): country, haul_id, ices_rect, month_year, substrate
dbl (21): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed 2 variables (X, Y)Make the grid with depth
First make a grid for the biomass data, then subset that based on the extend of the stomach data
x <- d$X
y <- d$Y
z <- chull(x, y)
coords <- cbind(x[z], y[z])
coords <- rbind(coords, coords[1, ])
plot(coords[, 1] ~ coords[, 2]) # plot data
sp_poly <- sp::SpatialPolygons(
list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
)
sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
data = data.frame(ID = 1)
)
cell_width <- 3
pred_grid <- expand.grid(
X = seq(min(d$X), max(d$X), cell_width),
Y = seq(min(d$Y), max(d$Y), cell_width),
year = c(1993:2021)
)
ggplot(pred_grid |> filter(year == 2019), aes(X, Y)) +
geom_point(size = 0.1) +
theme_void() +
coord_sf()filter: removed 882,420 rows (97%), 31,515 rows remaining
sp::coordinates(pred_grid) <- c("X", "Y")
inside <- !is.na(sp::over(pred_grid, as(sp_poly_df, "SpatialPolygons")))
pred_grid <- pred_grid[inside, ]
pred_grid <- as.data.frame(pred_grid)
ggplot(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000)) +
geom_point(size = 0.001, alpha = 0.5) +
NULLfilter: removed 531,776 rows (97%), 18,992 rows remaining
plot_map +
geom_point(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
NULLfilter: removed 531,776 rows (97%), 18,992 rows remainingWarning: Removed 704 rows containing missing values (`geom_point()`).
# Add lat and lon
# Need to go from UTM to lat long for this one...
# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
xy <- as.matrix(pred_grid |> dplyr::select(X, Y) |> mutate(X = X*1000, Y = Y*1000))mutate: changed 550,768 values (100%) of 'X' (0 new NA) changed 550,768 values (100%) of 'Y' (0 new NA)v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84 +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]
pred_grid$lon <- lonlat[, 1]
pred_grid$lat <- lonlat[, 2]
ggplot(filter(pred_grid, year == 1999), aes(lon, lat)) + geom_point()filter: removed 531,776 rows (97%), 18,992 rows remaining
# Add depth now to remove islands and remaining land
# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
# https://emodnet.ec.europa.eu/geoviewer/
dep_raster <- terra::rast(paste0(home, "/data/Mean depth natural colour (with land).nc"))
class(dep_raster)[1] "SpatRaster"
attr(,"package")
[1] "terra"crs(dep_raster, proj = TRUE)[1] "+proj=longlat +datum=WGS84 +no_defs"plot(dep_raster)
pred_grid$depth <- terra::extract(dep_raster, pred_grid |> dplyr::select(lon, lat))$elevation
ggplot(pred_grid, aes(lon, lat, color = depth*-1)) +
geom_point()
pred_grid$depth <- pred_grid$depth*-1
pred_grid <- pred_grid |> drop_na(depth)drop_na: removed 76,995 rows (14%), 473,773 rows remainingpred_grid |>
filter(year == 1999) |>
drop_na(depth) |>
#mutate(water = ifelse(depth < 0.00000001, "N", "Y")) |>
ggplot(aes(X*1000, Y*1000, fill = depth)) +
geom_raster() +
NULLfilter: removed 457,436 rows (97%), 16,337 rows remainingdrop_na: no rows removed
plot_map +
geom_point(data = pred_grid, aes(X*1000, Y*1000), size = 0.001) +
geom_sf()Warning: Removed 20329 rows containing missing values (`geom_point()`).
plot_map +
geom_raster(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000, fill = depth), size = 0.001) +
geom_sf()filter: removed 457,436 rows (97%), 16,337 rows remainingWarning in geom_raster(data = filter(pred_grid, year == 1999), aes(X * 1000, :
Ignoring unknown parameters: `size`Warning: Removed 762 rows containing missing values (`geom_raster()`).
Substrate
First add lat and lon based on X and Y (utm)
substrate <- terra::rast(paste0(home, "/data/substrate_tif/BALANCE_SEABED_SEDIMENT.tif"))
newcrs <- "+proj=longlat +datum=WGS84"
substrate_longlat = terra::project(substrate, newcrs)
|---------|---------|---------|---------|
=========================================
plot(substrate_longlat)
# Now extract the values from the saduria raster to pred_grid
pred_grid$substrate <- terra::extract(substrate_longlat, pred_grid |> dplyr::select(lon, lat))$BALANCE_SEABED_SEDIMENT
unique(pred_grid$substrate) [1] 3.000000 2.207703 3.000000 2.000000 5.000000 2.751267 2.411566 3.000000
[9] NA 3.023438 3.832672 4.693979 4.816533 2.407200 2.675354 5.000000
[17] 2.515299 2.104980 5.000000 3.309301 2.215659 4.000000 2.418371 3.475651
[25] 3.136719 3.112971 2.414551 4.865684 2.082968 3.555176 3.868312 3.925187
[33] 3.036916 4.963260 2.703572 4.166504 3.295410 3.913574 2.989511 2.172908
[41] 3.014659 3.177487 4.189453 4.184814 2.021204 3.047852 2.189942 4.947398
[49] 2.084473 2.895818 4.771301 3.115845 2.466890 4.691040 2.189209 2.446289
[57] 3.239676 3.565918 2.035608 4.346191 4.295410 2.221002 3.051758 2.749512
[65] 3.481570 2.114804 2.060020 4.481906 3.367188 3.831543 1.000000 4.013517
[73] 3.182683 2.004699 4.814331 2.606201 2.217144 1.983107 3.768274 2.091238
[81] 3.106734 2.150879 3.954897 2.124500 3.680969 3.013173 3.387896 2.451569
[89] 2.987671 3.004656 1.230957 3.087122 4.163172 4.331752 2.847656 3.770158
[97] 3.125428 4.137214 3.885498 3.675119 4.481557 3.618304 3.535532 4.465820
[105] 2.394043 3.985858 3.518455 2.930088 4.570079 3.035331 4.109087 3.688005
[113] 3.070784 3.621094 2.227002 3.092285 4.952732 4.725516 3.052990 3.999512
[121] 4.006014 3.395858 3.495117 3.194824 2.990967 4.053308 3.811279 4.897461
[129] 4.969238 3.747452 4.660156 4.250894 2.874023 2.834995 2.995728 2.603516
[137] 4.352661 4.529053 2.199062 3.234210 4.856023 4.949766 3.394751 3.105818
[145] 4.118652 2.052490 3.812706 2.458496 3.877447 4.508667 4.997460 4.851562
[153] 3.038818 2.274655 2.144043 2.438636 3.721595 3.773949 3.876709 4.546921
[161] 4.142808 3.771973 3.371582 4.526794 4.744454 2.290039 2.186924 3.043457
[169] 3.958649 4.084717 2.060291 4.713074 4.014312 2.630943 4.024330 4.889182
[177] 4.974456 3.012340 4.260148 4.845399 3.074554 4.409346 3.289751 3.797133
[185] 4.741160 4.990303 4.992432 2.236874 4.020061 4.972815 3.099927 3.104312
[193] 3.131327 4.014556 4.802023 3.626754 2.661095 4.347060 3.388261 2.413086
[201] 3.486816 4.048115 4.273254 3.218506 4.880537 2.026447 3.734767 3.617798
[209] 3.750732 3.903769 3.027108 4.973728 2.000000 3.153192 4.907196 4.142573
[217] 2.127777 3.366211 3.833305 3.686989 3.969047 2.451714 3.201564 3.003133
[225] 3.956664 2.793945 3.498962 3.130957 3.051147 3.170654 3.023269 2.039551
[233] 4.997359 3.213900 2.033055 2.056885 3.334961 2.716752 2.680624 3.365484
[241] 4.719290 4.165028 3.439011 2.860311 4.233765 2.220473 3.163818 4.547607
[249] 2.650146 3.881714 2.063829 2.739014 2.771362 3.844120 4.561007 2.790447
[257] 3.166412 2.996099 3.074219 4.000000 4.284601 3.897705 2.213754 2.750740
[265] 3.372153 3.138967 4.034563 2.040667 3.311523 4.179138 2.755887 2.779238
[273] 3.002139 2.831177 4.645216 4.891490 4.455050 3.114136 2.985913 4.142618
[281] 4.969971 4.160156 2.700928 2.437500 3.261470 2.720773 2.674125 4.945322
[289] 4.977282 4.981665 2.354759 3.050781 3.362305 4.034168 3.144043 3.427612
[297] 3.116003 4.706665 4.172024 4.919174 4.779785 4.999612 3.258072 4.724852
[305] 2.403809 2.568726 4.936097 4.244507 4.219971 4.019287 4.826586 3.978870
[313] 4.504028 3.343700 4.380371 2.978496 4.912631 3.740993 3.936610 2.325874
[321] 4.247256 2.162385 2.233503 2.387978 4.661255 3.199463 3.273904 2.430664
[329] 2.529750 2.882263 3.927612 4.913956 2.105009 3.039979 3.861948 4.729696
[337] 3.608606 2.118617 2.201874 3.994564 2.977833 2.230173 3.841385 2.345091
[345] 1.502542 3.117752 3.326479 4.573112 2.669289 4.844727 3.000776 2.726397
[353] 3.340427 4.578949 4.333318 4.793458 1.626249 2.945679 3.053248 3.968295
[361] 2.776123 1.726807 3.005692 3.789629 2.304810 3.019462 4.060608 2.094116
[369] 1.936951 3.144531 3.184711 2.283091 2.409424 1.728929 4.570679 4.402838
[377] 4.916132 2.350037 4.675385 4.664328 3.292003 2.157376 3.532034 4.425058
[385] 4.074768 4.896721 4.358123 2.311279 3.726698 3.015447 3.681493 4.969622
[393] 2.488892 2.328379 3.550172 4.602893 3.802612 4.695945 4.574463 4.006090
[401] 3.893816 4.331733 3.730713 4.537598 2.354067 4.375757 4.304613 4.836199
[409] 4.999268 4.829123 3.318359 4.121094 2.585449 4.685639 4.022028 4.823242
[417] 4.928812 2.574462 2.792725 4.440335 4.267953 4.459229 4.037476 2.681763
[425] 3.636063 4.099476 3.484864 4.835382 4.030147 4.965707 3.045144 4.963444
[433] 4.151367 4.356689 4.242188 4.018444 4.086391 4.748169 4.875027 2.393384
[441] 4.611477 3.831052 3.592084 2.666593 4.245469 4.447388 3.375265 4.155447
[449] 3.733586 2.018555 4.718811 4.296753 3.305542 4.881466 4.944336 4.507011
[457] 4.219604 4.868286 2.236110 4.212211 4.623746 3.489990 4.505832factor(sort(unique(round(pred_grid$substrate))))[1] 1 2 3 4 5
Levels: 1 2 3 4 5pred_grid$substrate <- round(pred_grid$substrate)
pred_grid <- pred_grid %>% mutate(substrate = ifelse(substrate == 1, "bedrock", substrate),
substrate = ifelse(substrate == 2, "hard-bottom complex", substrate),
substrate = ifelse(substrate == 3, "sand", substrate),
substrate = ifelse(substrate == 4, "hard clay", substrate),
substrate = ifelse(substrate == 5, "mud", substrate))mutate: converted 'substrate' from double to character (0 new NA)# I. Bedrock.
# II. Hard bottom complex, includes patchy hard surfaces and coarse sand (sometimes also clay) to boulders.
# III. Sand including fine to coarse sand (with gravel exposures).
# IV. Hard clay sometimes/often/possibly exposed or covered with a thin layer of sand/gravel.
# V. Mud including gyttja-clay to gyttja-silt.
# Plot
ggplot(pred_grid, aes(X, Y, fill = substrate)) +
geom_raster() +
coord_sf()
Add month_year variable for matching with raster layers
pred_grid_1 <- pred_grid |> mutate(quarter = 1)mutate: new variable 'quarter' (double) with one unique value and 0% NApred_grid_4 <- pred_grid |> mutate(quarter = 4)mutate: new variable 'quarter' (double) with one unique value and 0% NAdat <- bind_rows(pred_grid_1, pred_grid_4) |>
mutate(month = ifelse(quarter == 1, 2, 11), # most common months
month_year = paste(month, year, sep = "_"))mutate: new variable 'month' (double) with 2 unique values and 0% NA new variable 'month_year' (character) with 58 unique values and 0% NAOxygen
# Downloaded from here: https://data.marine.copernicus.eu/products
# Extract raster points: https://gisday.wordpress.com/2014/03/24/extract-raster-values-from-points-using-r/comment-page-1/
# https://rpubs.com/boyerag/297592
# https://pjbartlein.github.io/REarthSysSci/netCDF.html#get-a-variable
# Open the netCDF file
ncin <- nc_open(paste0(home, "/data/NEMO/cmems_mod_bal_bgc_my_P1M-m_1691065496556.nc"))
print(ncin)File /Users/maxlindmark/Dropbox/Max work/R/spatial-metabolic-index/data/NEMO/cmems_mod_bal_bgc_my_P1M-m_1691065496556.nc (NC_FORMAT_CLASSIC):
1 variables (excluding dimension variables):
float o2b[lon,lat,time]
standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
long_name: Dissolved_Oxygen_Concentration
units: mmol m-3
_FillValue: -999
missing_value: -999
cell_methods: time: mean
unit_long: millimole O2 per unit volume
interval_write: 1 d
interval_operation: 90 s
online_operation: average
_ChunkSizes: 1
_ChunkSizes: 774
_ChunkSizes: 763
3 dimensions:
time Size:347
standard_name: time
long_name: time
bounds: time_bnds
units: days since 1900-01-01 00:00:00
calendar: standard
axis: T
_ChunkSizes: 512
_CoordinateAxisType: Time
valid_min: 34013
valid_max: 44544.5
lat Size:498
standard_name: latitude
long_name: Latitude of each location
units: degrees_north
axis: Y
_CoordinateAxisType: Lat
valid_min: 53.0082969665527
valid_max: 61.2914962768555
lon Size:519
standard_name: longitude
long_name: Longitude of each location
units: degrees_east
axis: X
_CoordinateAxisType: Lon
valid_min: 10.5414848327637
valid_max: 24.9298248291016
22 global attributes:
CDI: Climate Data Interface version 1.9.9rc1 (https://mpimet.mpg.de/cdi)
Conventions: CF-1.0
history: Tue Mar 14 13:31:41 2023: cdo -z zip_1 --sortname copy /net/isilon/ifs/arch/home/iri/BALMFC/Nemo4.0/BALMFCMYP2022_univariateSST//BAL-ERGOM_BGC-DailyMeans-20211201.nc tmp_20211201.nc
source: CMEMS BAL MFC NEMO model output converted to NetCDF
institution: Baltic MFC, PU Danish Meteorological Institute
comment: Data on cropped native product grid. Horizontal velocities destagged
grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
contact: servicedesk.cmems@mercator-ocean.eu
references: https://marine.copernicus.eu/
file_quality_index: 1
compression: yes
stop_date: 2021-12-01 12:00:00
creation_date: 2023-03-14 13:31:27
FROM_ORIGINAL_FILE__netcdf_version_id: 4.6.1 of Oct 19 2018 20:03:53 $
FROM_ORIGINAL_FILE__easternmost_longitude: 30.2074012756348
FROM_ORIGINAL_FILE__northernmost_latitude: 65.8909912109375
FROM_ORIGINAL_FILE__westernmost_longitude: 9.04154205322266
FROM_ORIGINAL_FILE__southernmost_latitude: 53.0082969665527
start_date: 2021-12-01 12:00:00
CDO: Climate Data Operators version 1.9.9rc1 (https://mpimet.mpg.de/cdo)
title: CMEMS ERGOM monthly integrated model fields
_CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)[1] 10.54148 10.56926 10.59704 10.62481 10.65259 10.68037lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163# Get time
time <- ncvar_get(ncin,"time")
time [1] 34013.0 34042.5 34073.0 34103.5 34134.0 34164.5 34195.5 34226.0 34256.5
[10] 34287.0 34317.5 34348.5 34378.0 34407.5 34438.0 34468.5 34499.0 34529.5
[19] 34560.5 34591.0 34621.5 34652.0 34682.5 34713.5 34743.0 34772.5 34803.0
[28] 34833.5 34864.0 34894.5 34925.5 34956.0 34986.5 35017.0 35047.5 35078.5
[37] 35108.5 35138.5 35169.0 35199.5 35230.0 35260.5 35291.5 35322.0 35352.5
[46] 35383.0 35413.5 35444.5 35474.0 35503.5 35534.0 35564.5 35595.0 35625.5
[55] 35656.5 35687.0 35717.5 35748.0 35778.5 35809.5 35839.0 35868.5 35899.0
[64] 35929.5 35960.0 35990.5 36021.5 36052.0 36082.5 36113.0 36143.5 36174.5
[73] 36204.0 36233.5 36264.0 36294.5 36325.0 36355.5 36386.5 36417.0 36447.5
[82] 36478.0 36508.5 36539.5 36569.5 36599.5 36630.0 36660.5 36691.0 36721.5
[91] 36752.5 36783.0 36813.5 36844.0 36874.5 36905.5 36935.0 36964.5 36995.0
[100] 37025.5 37056.0 37086.5 37117.5 37148.0 37178.5 37209.0 37239.5 37270.5
[109] 37300.0 37329.5 37360.0 37390.5 37421.0 37451.5 37482.5 37513.0 37543.5
[118] 37574.0 37604.5 37635.5 37665.0 37694.5 37725.0 37755.5 37786.0 37816.5
[127] 37847.5 37878.0 37908.5 37939.0 37969.5 38000.5 38030.5 38060.5 38091.0
[136] 38121.5 38152.0 38182.5 38213.5 38244.0 38274.5 38305.0 38335.5 38366.5
[145] 38396.0 38425.5 38456.0 38486.5 38517.0 38547.5 38578.5 38609.0 38639.5
[154] 38670.0 38700.5 38731.5 38761.0 38790.5 38821.0 38851.5 38882.0 38912.5
[163] 38943.5 38974.0 39004.5 39035.0 39065.5 39096.5 39126.0 39155.5 39186.0
[172] 39216.5 39247.0 39277.5 39308.5 39339.0 39369.5 39400.0 39430.5 39461.5
[181] 39491.5 39521.5 39552.0 39582.5 39613.0 39643.5 39674.5 39705.0 39735.5
[190] 39766.0 39796.5 39827.5 39857.0 39886.5 39917.0 39947.5 39978.0 40008.5
[199] 40039.5 40070.0 40100.5 40131.0 40161.5 40192.5 40222.0 40251.5 40282.0
[208] 40312.5 40343.0 40373.5 40404.5 40435.0 40465.5 40496.0 40526.5 40557.5
[217] 40587.0 40616.5 40647.0 40677.5 40708.0 40738.5 40769.5 40800.0 40830.5
[226] 40861.0 40891.5 40922.5 40952.5 40982.5 41013.0 41043.5 41074.0 41104.5
[235] 41135.5 41166.0 41196.5 41227.0 41257.5 41288.5 41318.0 41347.5 41378.0
[244] 41408.5 41439.0 41469.5 41500.5 41531.0 41561.5 41592.0 41622.5 41653.5
[253] 41683.0 41712.5 41743.0 41773.5 41804.0 41834.5 41865.5 41896.0 41926.5
[262] 41957.0 41987.5 42018.5 42048.0 42077.5 42108.0 42138.5 42169.0 42199.5
[271] 42230.5 42261.0 42291.5 42322.0 42352.5 42383.5 42413.5 42443.5 42474.0
[280] 42504.5 42535.0 42565.5 42596.5 42627.0 42657.5 42688.0 42718.5 42749.5
[289] 42779.0 42808.5 42839.0 42869.5 42900.0 42930.5 42961.5 42992.0 43022.5
[298] 43053.0 43083.5 43114.5 43144.0 43173.5 43204.0 43234.5 43265.0 43295.5
[307] 43326.5 43357.0 43387.5 43418.0 43448.5 43479.5 43509.0 43538.5 43569.0
[316] 43599.5 43630.0 43660.5 43691.5 43722.0 43752.5 43783.0 43813.5 43844.5
[325] 43874.5 43904.5 43935.0 43965.5 43996.0 44026.5 44057.5 44088.0 44118.5
[334] 44149.0 44179.5 44210.5 44240.0 44269.5 44300.0 44330.5 44361.0 44391.5
[343] 44422.5 44453.0 44483.5 44514.0 44544.5tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt[1] 347tunits$hasatt
[1] TRUE
$value
[1] "days since 1900-01-01 00:00:00"# Get oxygen
dname <- "o2b"
oxy_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(oxy_array)[1] 519 498 347# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")
# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])
# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))
# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)
# Replace netCDF fill values with NA's
oxy_array[oxy_array == fillvalue$value] <- NA
dim(oxy_array)[1] 519 498 347str(dim(oxy_array)) int [1:3] 519 498 347# The third slot is the date index
# Loop through all "dates", put into a list
dlist <- list()
for(i in 1:length(months)) {
oxy_sub <- oxy_array[, , i]
dlist[[i]] <- oxy_sub
}
# Name the list
names(dlist) <- paste(months, years, sep = "_")
str(dlist)List of 347
$ 2_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
[list output truncated]# Create data holding object
oxy_data_list <- list()
# Loop through each month_year and extract raster values for the cpue data points
for(i in unique(dat$month_year)) {
# Set plot limits
ymin = 54; ymax = 58; xmin = 12; xmax = 22
# Subset a month-year combination
oxy_slice <- dlist[[i]]
# Create raster for that year (i)
r <- raster(t(oxy_slice), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
# Flip...
r <- flip(r, direction = 'y')
plot(r, main = paste(i), ylim = c(ymin, ymax), xlim = c(xmin, xmax))
# Filter the same year (i) in the cpue data and select only coordinates
d_slice <- dat %>% filter(month_year == i) %>% dplyr::select(lon, lat)
# Make into a SpatialPoints object
data_sp <- SpatialPoints(d_slice)
# Extract raster value (oxygen)
rasValue <- raster::extract(r, data_sp)
# Now we want to plot the results of the raster extractions by plotting the cpue data points over a raster and saving it for each year.
# Make the SpatialPoints object into a raster again (for plot)
df <- as.data.frame(data_sp)
# Add in the raster value in the df holding the coordinates for the cpue data
d_slice$oxy <- rasValue
# Add in which year
d_slice$month_year <- i
# Now the unit of oxygen is mmol/m3. I want it to be ml/L. The original model is in unit ml/L
# and it's been converted by the data host. Since it was converted without accounting for
# pressure or temperature, I can simply use the following conversion factor:
# 1 ml/l = 103/22.391 = 44.661 μmol/l -> 1 ml/l = 0.044661 mmol/l = 44.661 mmol/m^3 -> 0.0223909 ml/l = 1mmol/m^3
# https://ocean.ices.dk/tools/unitconversion.aspx
d_slice$oxy <- d_slice$oxy * 0.0223909
# Add each years' data in the list
oxy_data_list[[i]] <- d_slice
}filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
# Now create a data frame from the list of all annual values
big_dat_oxy <- dplyr::bind_rows(oxy_data_list)Temperature
# Open the netCDF file
ncin <- nc_open(paste0(home, "/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065935335.nc"))
print(ncin)File /Users/maxlindmark/Dropbox/Max work/R/spatial-metabolic-index/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065935335.nc (NC_FORMAT_CLASSIC):
1 variables (excluding dimension variables):
float bottomT[lon,lat,time]
standard_name: sea_water_potential_temperature_at_sea_floor
long_name: Sea floor potential temperature
units: degrees_C
_FillValue: -999
missing_value: -999
cell_methods: time: mean
unit_long: degrees Celsius
interval_write: 1 d
interval_operation: 90 s
online_operation: average
_ChunkSizes: 1
_ChunkSizes: 774
_ChunkSizes: 763
3 dimensions:
time Size:347
standard_name: time
long_name: time
bounds: time_bnds
units: days since 1900-01-01 00:00:00
calendar: standard
axis: T
_ChunkSizes: 512
_CoordinateAxisType: Time
valid_min: 34013
valid_max: 44544.5
lat Size:498
standard_name: latitude
long_name: Latitude of each location
units: degrees_north
axis: Y
_CoordinateAxisType: Lat
valid_min: 53.0082969665527
valid_max: 61.2914962768555
lon Size:519
standard_name: longitude
long_name: Longitude of each location
units: degrees_east
axis: X
_CoordinateAxisType: Lon
valid_min: 10.5414848327637
valid_max: 24.9298248291016
22 global attributes:
CDI: Climate Data Interface version 1.9.9rc1 (https://mpimet.mpg.de/cdi)
Conventions: CF-1.0
source: CMEMS BAL MFC NEMO model output converted to NetCDF
institution: Baltic MFC, PU Danish Meteorological Institute
comment: Data on cropped native product grid. Horizontal velocities destagged
grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
contact: servicedesk.cmems@mercator-ocean.eu
references: https://marine.copernicus.eu/
file_quality_index: 1
compression: yes
stop_date: 2021-12-01 12:00:00
creation_date: 2023-03-14 13:30:28
FROM_ORIGINAL_FILE__netcdf_version_id: 4.6.1 of Oct 19 2018 20:03:53 $
FROM_ORIGINAL_FILE__easternmost_longitude: 30.2074012756348
FROM_ORIGINAL_FILE__northernmost_latitude: 65.8909912109375
FROM_ORIGINAL_FILE__westernmost_longitude: 9.04154205322266
FROM_ORIGINAL_FILE__southernmost_latitude: 53.0082969665527
start_date: 2021-12-01 12:00:00
CDO: Climate Data Operators version 1.9.9rc1 (https://mpimet.mpg.de/cdo)
title: CMEMS NEMO monthly integrated model fields
_CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
history: Data extracted from dataset http://localhost:8080/thredds/dodsC/cmems_mod_bal_phy_my_P1M-m# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)[1] 10.54148 10.56926 10.59704 10.62481 10.65259 10.68037lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163# Get time
time <- ncvar_get(ncin,"time")
time [1] 34013.0 34042.5 34073.0 34103.5 34134.0 34164.5 34195.5 34226.0 34256.5
[10] 34287.0 34317.5 34348.5 34378.0 34407.5 34438.0 34468.5 34499.0 34529.5
[19] 34560.5 34591.0 34621.5 34652.0 34682.5 34713.5 34743.0 34772.5 34803.0
[28] 34833.5 34864.0 34894.5 34925.5 34956.0 34986.5 35017.0 35047.5 35078.5
[37] 35108.5 35138.5 35169.0 35199.5 35230.0 35260.5 35291.5 35322.0 35352.5
[46] 35383.0 35413.5 35444.5 35474.0 35503.5 35534.0 35564.5 35595.0 35625.5
[55] 35656.5 35687.0 35717.5 35748.0 35778.5 35809.5 35839.0 35868.5 35899.0
[64] 35929.5 35960.0 35990.5 36021.5 36052.0 36082.5 36113.0 36143.5 36174.5
[73] 36204.0 36233.5 36264.0 36294.5 36325.0 36355.5 36386.5 36417.0 36447.5
[82] 36478.0 36508.5 36539.5 36569.5 36599.5 36630.0 36660.5 36691.0 36721.5
[91] 36752.5 36783.0 36813.5 36844.0 36874.5 36905.5 36935.0 36964.5 36995.0
[100] 37025.5 37056.0 37086.5 37117.5 37148.0 37178.5 37209.0 37239.5 37270.5
[109] 37300.0 37329.5 37360.0 37390.5 37421.0 37451.5 37482.5 37513.0 37543.5
[118] 37574.0 37604.5 37635.5 37665.0 37694.5 37725.0 37755.5 37786.0 37816.5
[127] 37847.5 37878.0 37908.5 37939.0 37969.5 38000.5 38030.5 38060.5 38091.0
[136] 38121.5 38152.0 38182.5 38213.5 38244.0 38274.5 38305.0 38335.5 38366.5
[145] 38396.0 38425.5 38456.0 38486.5 38517.0 38547.5 38578.5 38609.0 38639.5
[154] 38670.0 38700.5 38731.5 38761.0 38790.5 38821.0 38851.5 38882.0 38912.5
[163] 38943.5 38974.0 39004.5 39035.0 39065.5 39096.5 39126.0 39155.5 39186.0
[172] 39216.5 39247.0 39277.5 39308.5 39339.0 39369.5 39400.0 39430.5 39461.5
[181] 39491.5 39521.5 39552.0 39582.5 39613.0 39643.5 39674.5 39705.0 39735.5
[190] 39766.0 39796.5 39827.5 39857.0 39886.5 39917.0 39947.5 39978.0 40008.5
[199] 40039.5 40070.0 40100.5 40131.0 40161.5 40192.5 40222.0 40251.5 40282.0
[208] 40312.5 40343.0 40373.5 40404.5 40435.0 40465.5 40496.0 40526.5 40557.5
[217] 40587.0 40616.5 40647.0 40677.5 40708.0 40738.5 40769.5 40800.0 40830.5
[226] 40861.0 40891.5 40922.5 40952.5 40982.5 41013.0 41043.5 41074.0 41104.5
[235] 41135.5 41166.0 41196.5 41227.0 41257.5 41288.5 41318.0 41347.5 41378.0
[244] 41408.5 41439.0 41469.5 41500.5 41531.0 41561.5 41592.0 41622.5 41653.5
[253] 41683.0 41712.5 41743.0 41773.5 41804.0 41834.5 41865.5 41896.0 41926.5
[262] 41957.0 41987.5 42018.5 42048.0 42077.5 42108.0 42138.5 42169.0 42199.5
[271] 42230.5 42261.0 42291.5 42322.0 42352.5 42383.5 42413.5 42443.5 42474.0
[280] 42504.5 42535.0 42565.5 42596.5 42627.0 42657.5 42688.0 42718.5 42749.5
[289] 42779.0 42808.5 42839.0 42869.5 42900.0 42930.5 42961.5 42992.0 43022.5
[298] 43053.0 43083.5 43114.5 43144.0 43173.5 43204.0 43234.5 43265.0 43295.5
[307] 43326.5 43357.0 43387.5 43418.0 43448.5 43479.5 43509.0 43538.5 43569.0
[316] 43599.5 43630.0 43660.5 43691.5 43722.0 43752.5 43783.0 43813.5 43844.5
[325] 43874.5 43904.5 43935.0 43965.5 43996.0 44026.5 44057.5 44088.0 44118.5
[334] 44149.0 44179.5 44210.5 44240.0 44269.5 44300.0 44330.5 44361.0 44391.5
[343] 44422.5 44453.0 44483.5 44514.0 44544.5tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt[1] 347tunits$hasatt
[1] TRUE
$value
[1] "days since 1900-01-01 00:00:00"# Get temperature
dname <- "bottomT"
temp_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(temp_array)[1] 519 498 347# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")
# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])
# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))
# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)
# Replace netCDF fill values with NA's
temp_array[temp_array == fillvalue$value] <- NA
# Loop through all "dates", put into a list
dlist <- list()
for(i in 1:length(months)) {
temp_sub <- temp_array[, , i]
dlist[[i]] <- temp_sub
}
# Name the list
names(dlist) <- paste(months, years, sep = "_")
str(dlist)List of 347
$ 2_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
[list output truncated]# Create data holding object
temp_data_list <- list()
# Loop through each month_year and extract raster values for the cpue data points
for(i in unique(dat$month_year)) { # We can use q1 as looping index, doesn't matter!
# Set plot limits
ymin = 54; ymax = 58; xmin = 12; xmax = 22
# Subset a month-year combination
temp_slice <- dlist[[i]]
# Create raster for that year (i)
r <- raster(t(temp_slice), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
# Flip...
r <- flip(r, direction = 'y')
plot(r, main = paste(i), ylim = c(ymin, ymax), xlim = c(xmin, xmax))
# Filter the same year (i) in the cpue data and select only coordinates
d_slice <- dat %>% filter(month_year == i) %>% dplyr::select(lon, lat)
# Make into a SpatialPoints object
data_sp <- SpatialPoints(d_slice)
# Extract raster value (oxygen)
rasValue <- raster::extract(r, data_sp)
# Now we want to plot the results of the raster extractions by plotting the cpue data points over a raster and saving it for each year.
# Make the SpatialPoints object into a raster again (for plot)
df <- as.data.frame(data_sp)
# Add in the raster value in the df holding the coordinates for the cpue data
d_slice$temp <- rasValue
# Add in which year
d_slice$month_year <- i
# Add each years' data in the list
temp_data_list[[i]] <- d_slice
}filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
# Now create a data frame from the list of all annual values
big_dat_temp <- dplyr::bind_rows(temp_data_list)Salinity
# Open the netCDF file
ncin <- nc_open(paste0(home, "/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065974147.nc"))
print(ncin)File /Users/maxlindmark/Dropbox/Max work/R/spatial-metabolic-index/data/NEMO/cmems_mod_bal_phy_my_P1M-m_1691065974147.nc (NC_FORMAT_CLASSIC):
1 variables (excluding dimension variables):
float sob[lon,lat,time]
standard_name: sea_water_salinity_at_sea_floor
long_name: Sea water salinity at sea floor
units: 0.001
_FillValue: -999
missing_value: -999
cell_methods: time: mean
unit_long: 0.001
interval_write: 1 d
interval_operation: 90 s
online_operation: average
_ChunkSizes: 1
_ChunkSizes: 774
_ChunkSizes: 763
3 dimensions:
time Size:347
standard_name: time
long_name: time
bounds: time_bnds
units: days since 1900-01-01 00:00:00
calendar: standard
axis: T
_ChunkSizes: 512
_CoordinateAxisType: Time
valid_min: 34013
valid_max: 44544.5
lat Size:498
standard_name: latitude
long_name: Latitude of each location
units: degrees_north
axis: Y
_CoordinateAxisType: Lat
valid_min: 53.0082969665527
valid_max: 61.2914962768555
lon Size:519
standard_name: longitude
long_name: Longitude of each location
units: degrees_east
axis: X
_CoordinateAxisType: Lon
valid_min: 10.5414848327637
valid_max: 24.9298248291016
22 global attributes:
CDI: Climate Data Interface version 1.9.9rc1 (https://mpimet.mpg.de/cdi)
Conventions: CF-1.0
source: CMEMS BAL MFC NEMO model output converted to NetCDF
institution: Baltic MFC, PU Danish Meteorological Institute
comment: Data on cropped native product grid. Horizontal velocities destagged
grid_resolution: ~1 nautical mile (1min latitude; 1min40sec longitude)
contact: servicedesk.cmems@mercator-ocean.eu
references: https://marine.copernicus.eu/
file_quality_index: 1
compression: yes
stop_date: 2021-12-01 12:00:00
creation_date: 2023-03-14 13:30:28
FROM_ORIGINAL_FILE__netcdf_version_id: 4.6.1 of Oct 19 2018 20:03:53 $
FROM_ORIGINAL_FILE__easternmost_longitude: 30.2074012756348
FROM_ORIGINAL_FILE__northernmost_latitude: 65.8909912109375
FROM_ORIGINAL_FILE__westernmost_longitude: 9.04154205322266
FROM_ORIGINAL_FILE__southernmost_latitude: 53.0082969665527
start_date: 2021-12-01 12:00:00
CDO: Climate Data Operators version 1.9.9rc1 (https://mpimet.mpg.de/cdo)
title: CMEMS NEMO monthly integrated model fields
_CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
history: Data extracted from dataset http://localhost:8080/thredds/dodsC/cmems_mod_bal_phy_my_P1M-m# Get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)[1] 10.54148 10.56926 10.59704 10.62481 10.65259 10.68037lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)[1] 53.00830 53.02496 53.04163 53.05830 53.07496 53.09163# Get time
time <- ncvar_get(ncin,"time")
time [1] 34013.0 34042.5 34073.0 34103.5 34134.0 34164.5 34195.5 34226.0 34256.5
[10] 34287.0 34317.5 34348.5 34378.0 34407.5 34438.0 34468.5 34499.0 34529.5
[19] 34560.5 34591.0 34621.5 34652.0 34682.5 34713.5 34743.0 34772.5 34803.0
[28] 34833.5 34864.0 34894.5 34925.5 34956.0 34986.5 35017.0 35047.5 35078.5
[37] 35108.5 35138.5 35169.0 35199.5 35230.0 35260.5 35291.5 35322.0 35352.5
[46] 35383.0 35413.5 35444.5 35474.0 35503.5 35534.0 35564.5 35595.0 35625.5
[55] 35656.5 35687.0 35717.5 35748.0 35778.5 35809.5 35839.0 35868.5 35899.0
[64] 35929.5 35960.0 35990.5 36021.5 36052.0 36082.5 36113.0 36143.5 36174.5
[73] 36204.0 36233.5 36264.0 36294.5 36325.0 36355.5 36386.5 36417.0 36447.5
[82] 36478.0 36508.5 36539.5 36569.5 36599.5 36630.0 36660.5 36691.0 36721.5
[91] 36752.5 36783.0 36813.5 36844.0 36874.5 36905.5 36935.0 36964.5 36995.0
[100] 37025.5 37056.0 37086.5 37117.5 37148.0 37178.5 37209.0 37239.5 37270.5
[109] 37300.0 37329.5 37360.0 37390.5 37421.0 37451.5 37482.5 37513.0 37543.5
[118] 37574.0 37604.5 37635.5 37665.0 37694.5 37725.0 37755.5 37786.0 37816.5
[127] 37847.5 37878.0 37908.5 37939.0 37969.5 38000.5 38030.5 38060.5 38091.0
[136] 38121.5 38152.0 38182.5 38213.5 38244.0 38274.5 38305.0 38335.5 38366.5
[145] 38396.0 38425.5 38456.0 38486.5 38517.0 38547.5 38578.5 38609.0 38639.5
[154] 38670.0 38700.5 38731.5 38761.0 38790.5 38821.0 38851.5 38882.0 38912.5
[163] 38943.5 38974.0 39004.5 39035.0 39065.5 39096.5 39126.0 39155.5 39186.0
[172] 39216.5 39247.0 39277.5 39308.5 39339.0 39369.5 39400.0 39430.5 39461.5
[181] 39491.5 39521.5 39552.0 39582.5 39613.0 39643.5 39674.5 39705.0 39735.5
[190] 39766.0 39796.5 39827.5 39857.0 39886.5 39917.0 39947.5 39978.0 40008.5
[199] 40039.5 40070.0 40100.5 40131.0 40161.5 40192.5 40222.0 40251.5 40282.0
[208] 40312.5 40343.0 40373.5 40404.5 40435.0 40465.5 40496.0 40526.5 40557.5
[217] 40587.0 40616.5 40647.0 40677.5 40708.0 40738.5 40769.5 40800.0 40830.5
[226] 40861.0 40891.5 40922.5 40952.5 40982.5 41013.0 41043.5 41074.0 41104.5
[235] 41135.5 41166.0 41196.5 41227.0 41257.5 41288.5 41318.0 41347.5 41378.0
[244] 41408.5 41439.0 41469.5 41500.5 41531.0 41561.5 41592.0 41622.5 41653.5
[253] 41683.0 41712.5 41743.0 41773.5 41804.0 41834.5 41865.5 41896.0 41926.5
[262] 41957.0 41987.5 42018.5 42048.0 42077.5 42108.0 42138.5 42169.0 42199.5
[271] 42230.5 42261.0 42291.5 42322.0 42352.5 42383.5 42413.5 42443.5 42474.0
[280] 42504.5 42535.0 42565.5 42596.5 42627.0 42657.5 42688.0 42718.5 42749.5
[289] 42779.0 42808.5 42839.0 42869.5 42900.0 42930.5 42961.5 42992.0 43022.5
[298] 43053.0 43083.5 43114.5 43144.0 43173.5 43204.0 43234.5 43265.0 43295.5
[307] 43326.5 43357.0 43387.5 43418.0 43448.5 43479.5 43509.0 43538.5 43569.0
[316] 43599.5 43630.0 43660.5 43691.5 43722.0 43752.5 43783.0 43813.5 43844.5
[325] 43874.5 43904.5 43935.0 43965.5 43996.0 44026.5 44057.5 44088.0 44118.5
[334] 44149.0 44179.5 44210.5 44240.0 44269.5 44300.0 44330.5 44361.0 44391.5
[343] 44422.5 44453.0 44483.5 44514.0 44544.5tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt[1] 347tunits$hasatt
[1] TRUE
$value
[1] "days since 1900-01-01 00:00:00"# Get Salinity
dname <- "sob"
sal_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(sal_array)[1] 519 498 347# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")
# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])
# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))
# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)
# Replace netCDF fill values with NA's
sal_array[sal_array == fillvalue$value] <- NA
# Loop through all "dates", put into a list
dlist <- list()
for(i in 1:length(months)) {
sal_sub <- sal_array[, , i]
dlist[[i]] <- sal_sub
}
# Name the list
names(dlist) <- paste(months, years, sep = "_")
str(dlist)List of 347
$ 2_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1993 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1993: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1994 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1994: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1995 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1995: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1996 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1996: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1997 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1997: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1998 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1998: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_1999 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_1999: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 5_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 6_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 7_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 8_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 9_2000 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 10_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 11_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 12_2000: num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 1_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 2_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 3_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
$ 4_2001 : num [1:519, 1:498] NA NA NA NA NA NA NA NA NA NA ...
[list output truncated]# Create data holding object
sal_data_list <- list()
# Loop through each month_year and extract raster values for the cpue data points
for(i in unique(dat$month_year)) { # We can use q1 as looping index, doesn't matter!
# Set plot limits
ymin = 54; ymax = 58; xmin = 12; xmax = 22
# Subset a month-year combination
sal_slice <- dlist[[i]]
# Create raster for that year (i)
r <- raster(t(sal_slice), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
# Flip...
r <- flip(r, direction = 'y')
plot(r, main = paste(i), ylim = c(ymin, ymax), xlim = c(xmin, xmax))
# Filter the same year (i) in the cpue data and select only coordinates
d_slice <- dat %>% filter(month_year == i) %>% dplyr::select(lon, lat)
# Make into a SpatialPoints object
data_sp <- SpatialPoints(d_slice)
# Extract raster value (oxygen)
rasValue <- raster::extract(r, data_sp)
# Now we want to plot the results of the raster extractions by plotting the cpue data points over a raster and saving it for each year.
# Make the SpatialPoints object into a raster again (for plot)
df <- as.data.frame(data_sp)
# Add in the raster value in the df holding the coordinates for the cpue data
d_slice$sal <- rasValue
# Add in which year
d_slice$month_year <- i
# Add each years' data in the list
sal_data_list[[i]] <- d_slice
}filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
filter: removed 931,209 rows (98%), 16,337 rows remaining
# Now create a data frame from the list of all annual values
big_dat_sal <- dplyr::bind_rows(sal_data_list)env_dat <- left_join(big_dat_oxy, big_dat_temp, by = c("month_year", "lon", "lat")) |>
left_join(big_dat_sal, by = c("month_year", "lon", "lat"))left_join: added one column (temp) > rows only in x 0 > rows only in y ( 0) > matched rows 947,546 > ========= > rows total 947,546left_join: added one column (sal) > rows only in x 0 > rows only in y ( 0) > matched rows 947,546 > ========= > rows total 947,546# Now join these data with the full_dat
dat_full <- left_join(dat, env_dat, by = c("month_year", "lon", "lat"))left_join: added 3 columns (oxy, temp, sal) > rows only in x 0 > rows only in y ( 0) > matched rows 947,546 > ========= > rows total 947,546ggplot(dat_full |> filter(quarter == 1), aes(X, Y, color = oxy)) +
geom_point() +
facet_wrap(~year)filter: removed 473,773 rows (50%), 473,773 rows remaining
Add ICES areas
# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile(paste0(home, "/data/ICES_StatRec_mapto_ICES_Areas/StatRec_map_Areas_Full_20170124.shp"))
head(shape) OBJECTID ID ICESNAME SOUTH WEST NORTH EAST AREA_KM2 stat_x stat_y Area_27
1 1 47 47A0 59.0 -44 59.5 -43 3178 -43.5 59.25 14.b.2
2 2 48 48A0 59.5 -44 60.0 -43 3132 -43.5 59.75 14.b.2
3 3 49 49A0 60.0 -44 60.5 -43 3085 -43.5 60.25 14.b.2
4 4 50 50A0 60.5 -44 61.0 -43 3038 -43.5 60.75 14.b.2
5 5 51 51A0 61.0 -44 61.5 -43 2991 -43.5 61.25 14.b.2
6 6 52 52A0 61.5 -44 62.0 -43 2944 -43.5 61.75 14.b.2
Perc MaxPer RNDMaxPer AreasList Shape_Leng Shape_Area
1 100.00000000 100.0000000 100 14.b.2 3 0.5
2 84.12674145 84.1267414 84 14.b.2 3 0.5
3 24.99803694 24.9980369 25 14.b.2 3 0.5
4 11.97744244 11.9774424 12 14.b.2 3 0.5
5 3.89717625 3.8971762 4 14.b.2 3 0.5
6 0.28578428 0.2857843 0 14.b.2 3 0.5pts <- SpatialPoints(cbind(dat_full$lon, dat_full$lat),
proj4string = CRS(proj4string(shape)))
dat_full$subdiv <- over(pts, shape)$Area_27
# Rename subdivisions to the more common names and do some more filtering (by sub div and area)
sort(unique(dat_full$subdiv))[1] "3.d.24" "3.d.25" "3.d.26" "3.d.27" "3.d.28.2"dat_full <- dat_full |>
mutate(sub_div = factor(subdiv),
sub_div = fct_recode(subdiv,
"24" = "3.d.24",
"25" = "3.d.25",
"26" = "3.d.26",
"27" = "3.d.27",
"28" = "3.d.28.1",
"28" = "3.d.28.2",
"29" = "3.d.29"),
sub_div = as.character(sub_div)) |>
filter(sub_div %in% c("24", "25", "26", "27", "28", 2)) |>
filter(lat > 54 & lat < 59 & lon < 22)Warning: There was 1 warning in `.fun()`.
ℹ In argument: `sub_div = fct_recode(...)`.
Caused by warning:
! Unknown levels in `f`: 3.d.28.1, 3.d.29mutate: new variable 'sub_div' (character) with 5 unique values and 0% NAfilter: no rows removed
filter: no rows removed# Add ICES rectangles
dat_full$ices_rect <- mapplots::ices.rect2(lon = dat_full$lon, lat = dat_full$lat)
plot_map +
geom_raster(data = filter(dat_full, year == 1999), aes(X*1000, Y*1000, fill = oxy)) +
facet_wrap(~sub_div)filter: removed 914,872 rows (97%), 32,674 rows remainingWarning: Removed 1524 rows containing missing values (`geom_raster()`).
dat_full <- dat_full |> dplyr::select(-subdiv)Save
# Remove variables and save
pred_grid_93_06 <- dat_full |> filter(year < 2007)filter: removed 490,110 rows (52%), 457,436 rows remainingpred_grid_07_19 <- dat_full |> filter(year > 2006)filter: removed 457,436 rows (48%), 490,110 rows remainingwrite_csv(pred_grid_93_06, file = paste0(home, "/data/clean/pred_grid_(1_2).csv"))
write_csv(pred_grid_07_19, file = paste0(home, "/data/clean/pred_grid_(2_2).csv"))knitr::knit_exit()